Non-Oscillatory  Central  Schemes  for  One-  and  Two-Dimensional 
MHD  Equations.  II:  High-order  Semi-discrete  Schemes* 

Jorge  Baibas'  and  Eitan  Tadmoh 
June  20,  2004 


Abstract 

We  present  a  new  family  of  high-resolution,  non-oscillatory  semi-discrete  central  schemes  for  the  approx¬ 
imate  solution  of  the  ideal  Magnetohydrodynamics  (MHD)  equations.  This  is  the  second  part  of  our  work, 
where  we  are  passing  from  the  fully-discrete,  staggered  schemes  utilized  in  [4]  to  a  semi-discrete  formulation 
advocated  in  [17].  This  semi-discrete  formulation  retains  the  simplicity  of  fully-discrete  central  schemes  while 
enhancing  efficency  and  adding  versatility.  The  semi-discrete  algorithm  offers  a  wider  range  of  options  to 
implement  its  two  key  steps:  non-oscillatory  reconstruction  and  evolution.  Along  with  the  description  of  the 
numerical  methods  employed,  we  present  several  prototype  MHD  problems.  Solutions  of  one-dimensional 
shock-tube  problems  and  the  two-dimensional  Kelvin-Helmholtz  and  Orszag-Tang  problems  are  carried  out 
using  third-  and  fourth-order  central  schemes  based  on  the  CWENO  reconstructions  of  Kurganov  and  Levy, 

[15],  and  Levy  et.  ah,  [21].  These  results  complement  those  presented  in  [4]  (and  the  references  therein)  and 
confirm  the  remarkable  versatility  of  central  schemes  as  black-box,  Jacobian-free  MHD  solvers. 

AMS  subject  classification:  Primary  65M10;  Secondary  65M05,  76W05 

Key  words.  Multidimensional  conservation  laws,  ideal  Magnetohydrodynamics  (MHD)  equations,  high-resolution 
central  schemes,  non-oscillatory  reconstructions,  Jacobian-free  form,  semi-discrete  schemes. 


1  Introduction 

In  this  paper  we  present  third-  and  fourth-order  accurate,  non-oscillatory  central  schemes  for  the  approximate 
solution  of  the  equations  of  ideal  Magnetohydrodynamics, 
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Here,  p  and  e  are  scalar  quantities  representing  respectively,  the  mass  density  and  the  total  internal  energy, 
v  =  ( vx,vy,vz)T  is  the  velocity  field  with  Euclidean  norm  v 2  :=  ||v||2,  and  B  =  (Bx,  By,  BZ)T  and  B2  :=  ||Bj|2 
represent  the  magnetic  field  and  its  norm.  The  pressure,  p,  is  coupled  to  the  internal  energy,  e  =  \pv2  +  ^B2  + 
p/(y  —  1),  where  7  is  the  (fixed)  ratio  of  specific  heats.  The  system  is  augmented  by  the  solenoidal  constraint, 

V  •  B  =  0;  (1.5) 

that  is,  if  the  condition  V  •  B  =  0  is  satisfied  initially  at  t  =  0,  then  by  (1.3)  it  remains  invariant  in  time. 
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Our  work  [4]  demonstrated  the  capability  of  central  schemes  to  accurately  and  efficiently  detect  and  resolve 
the  steep  gradients  that  characterize  the  solutions  of  (1.1)-(1.4);  it  suggested  further  development  of  higher  order, 
black-box  type  central  schemes  for  these  equations.  In  this  paper,  we  present  a  family  of  third-  and  fourth- 
order  accurate,  non-oscillatory  semi-discrete  central  schemes.  This  semi-discrete  formulation,  [17],  retains  the 
main  advantages  of  central  schemes  over  the  upwind  ones,  namely,  simplicity  and  efficiency  resulting  from 
using  a  minimal  amount  of  characteristic  information  when  evolving  the  solution  over  staggered  grids,  the 
use  of  Jacobian-free  formulations  and  avoiding  dimensional  splitting  in  multidimensional  models.  Moreover, 
for  r  -  order  methods,  the  numerical  dissipation  of  the  semi-discrete  formulation  is  of  order  0((Ax)2r~1)  vs. 
the  corresponding  fully-discrete  version  of  order  0((Ax)2r /At);  this  allows  the  use  of  smaller  time  steps  with 
considerably  less  smearing  of  discontinuities.  Finally,  the  semi-discrete  formulation  also  enjoys  a  remarkable 
adaptability  to  a  wide  range  of  non-oscillatory  reconstruction  techniques  and  evolution  routines  developed 
independently  of  the  methods. 

We  implement  these  new  schemes  to  approximate  the  solution  of  the  five  MHD  prototype  problems  discussed 
in  [4]:  two  one-dimensional  Brio  and  Wu  shock  tube  models  originally  proposed  in  [6],  two  different  configurations 

—  periodic  and  convective  models  of  the  2D  transverse  Kelvin-Helmholtz  instability  problem,  and  the  2D  MHD 
vortex  system  of  Orszag  and  Tang,  whose  evolution  involves  the  interaction  between  several  shock  waves  traveling 
at  various  speed  regimes,  making  it  a  very  suitable  model  to  test  the  robustness  of  numerical  algorithms. 

The  paper  is  structured  as  follows: 

In  §2,  we  extend  the  derivation  of  the  second-  and  third-order,  central  semi-discrete  schemes  presented 
in  [17,  16,  15]  to  schemes  of  arbitrary  order  r.  Starting  with  the  staggered  discretization  that  led  to  the 
fully-discrete  central  schemes,  [23,  22,  13,  20],  we  evolve  separately  the  smooth  and  non-smooth  parts  of  the 
solution  and  project  the  resulting  cell  averages  back  onto  the  original  grid.  We  then  pass  to  the  At  J,  0-limit, 
yielding  the  rth- order  semi-discrete  formulation  of  central  schemes.  This  rather  general  formulation  of  central 
schemes  entertains  a  broad  range  of  options  to  implement  the  two  main  steps  of  the  algorithm:  non-oscillatory 
reconstruction  of  point  values  from  cell  averages,  followed  by  the  evolution  of  the  corresponding  point  valued 
fluxes.  In  §3  we  discuss  the  non-oscillatory  reconstruction  techniques  and  evolution  routines  that  we  implement 
to  approximate  the  semi-discrete  formulation  of  equations  (1.1) — (1.4) .  For  the  reconstruction  of  point  values, 
we  propose  an  extension  of  the  third-order  Central  WENO  (CWENO)  reconstruction  of  Kurganov  and  Levy, 
[15],  and  the  genuinely  two-dimensional  fourth-order  CWENO  reconstruction  of  Levy  et.  al.,  [21],  developed 
originally  in  the  context  of  fully-discrete  central  schemes.  For  the  evolution  step,  we  use  third-  and  fourth-order 
strong  stability  preserving  (SSP)  Runge-Kutta  schemes,  [26,  28,  8].  We  complete  the  description  of  the  methods 
with  a  discussion  on  the  divergence-free  condition  of  the  magnetic  field,  V  •  B  =  0.  This  constraint  is  guaranteed 
in  the  one-dimensional  shock  tube  problems  and  the  two-dimensional  Kelvin-Helmholtz  instability  models,  but 
not  in  the  Orszag-Tang  MHD  vortex  system.  In  the  latter  case,  a  Leray  projection  corrector  is  implemented  at 
the  end  of  each  time-step,  replacing  the  computed  magnetic  field  with  its  divergence  free  projection. 

In  §4  and  §5  we  describe  the  prototype  MHD  problems  mentioned  above  and  present  the  numerical  solutions 
obtained  with  our  new  family  of  high  resolution  semi-discrete  central  schemes.  In  §4  we  report  on  the  one¬ 
dimensional  shock  tube  problems  of  Brio  and  Wu  [6]  and  discuss  the  improved  resolution  and  better  control 
of  spurious  oscillations  provided  by  the  new  schemes.  The  two-dimensional  models  and  the  corresponding 
numerical  results  are  presented  and  discussed  in  §5.  These  results  confirm  the  efficiency  of  central  schemes  to 
accurately  compute  the  solution  of  MHD  equations.  The  combination  of  simplicity  and  high  resolution  allows 
us  to  obtain  solutions  comparable  to  those  obtained  with  higher  order  upwind  schemes,  e.g.,  [6,  7,  14,  25,  24], 
without  having  to  refine  the  computational  mesh  as  in  [4] .  We  therefore  retain  the  simplicity  of  central  schemes 

—  they  eliminate  the  need  for  a  detailed  knowledge  of  the  eigen-structure  of  the  Jacobian  matrices  and  avoid 
dimensional  splitting  and  the  costly  use  of  (approximate)  Riemann  solvers  that  serve  as  building  blocks  for 
upwind  schemes;  at  the  same  time,  this  simplicity  is  achieved  without  an  increase  in  the  computational  work 
that  was  enforced  by  the  finer  meshes  required  with  the  lower  order  methods. 

2  Semi-discrete  Central  Schemes 

In  [4]  we  demonstrated  the  robustness  of  black-box  type  central  schemes  to  calculate  the  approximate  solutions 
of  equations  (1.1)- (1.4)  and  discussed  the  advantages  of  these  methods  over  those  based  on  upwind  differencing. 
We  also  pointed  out,  however,  that  the  efficiency  resulting  from  a  minimal  amount  of  characteristic  information 
required  by  central  schemes,  may  be  compromised  by  the  low-order  accuracy  of  the  schemes,  since  the  latter  often 
require  finer  meshes  to  resolve  steep  gradients.  In  order  to  circumvent  this  limitation  and  gain  full  advantage 
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of  the  simplicity,  ease  of  implementation  and  efficiency  of  central  schemes,  we  seek  higher  order  methods. 
Although  it  is  possible  to  carry  out  a  higher  order  extension  of  the  fully-discrete  central  schemes  discussed  in 
[4],  see  e.g.,  [19,  21],  we  find  the  semi-discrete  formulation  of  Kurganov  and  Tadmor,  [17],  more  advantageous: 
it  enjoys  a  reduced  numerical  dissipation  (of  order  0(( Ax)2r_1))  in  comparison  to  that  of  fully-discrete  schemes 
(~  0((Ax)2r  /  At)) ,  and  it  is  derived  independently  of  any  reconstruction  or  evolution  algorithm.  Indeed, 
many  of  the  existing  high-order,  non-oscillatory  reconstruction  algorithms  and  evolution  routines  can  be  easily 
incorporated  into  the  semi-discrete  central  formulation,  making  the  simplicity  and  ease  of  implementation  of 
central  schemes  an  attractive  alternative. 

2.1  One-dimensional  Schemes 

The  solutions  of  non-linear  hyperbolic  systems  of  conservation  laws, 


ut  +  f(u)x  =  0, 


(2.1) 


and  in  particular  those  of  equations  (1.1)-(1.4),  are  characterized  by  the  onset  and  propagation  of  discontinuities. 
The  higher  order  extensions  of  the  first-order  Lax-Friedrichs  central  scheme,  [18],  that  have  been  developed  over 
the  past  two  decades,  see  e.g.,  [23,  22,  13,  17],  provide  a  rather  general  approach  for  solving  these  type  of 
problems  with  only  a  minimal  amount  of  information  on  the  eigen  structure  of  the  Jacobian  matrix  of  f{u) 
(indeed,  Jacobian- free  formulations  of  central  schemes  completely  avoid  the  computation  of  |^).  These  schemes 

are  based  on  the  evolution  of  the  cell  averages  (we  use  the  usual  j-  to  denote  the  average,  1/|/|  •  J), 


u ‘  := 


:=  J- 1  u(£,tn)d£, 


T  Ax  Ax, 

d-a  ■ —  \p^a  5  J-  “ ^  J 


while  alternating  over  staggered  grids,  {Ij}  x  tn  and  {/J+i }  x  tn+l .  Dividing  (2.1)  by  Ax  and  integrating  over 
the  control  volume  x  [tn,tn+l]  yields  an  exact  discrete  statement  of  the  conservation  law  in  terms  of  these 
cell  averages, 


?"7ra+1  —  an 

Uj+ |  -  uo+h 


1 

Ax 


f(u(xj+i,T))dr  - 


f(u(xj,T))dr 


(2.2) 


The  evaluation  of  the  expressions  on  the  right  of  (2.2)  proceeds  in  two  steps. 

First,  the  cell  averages  {it"  1 }  are  recovered  by  integrating  a  piecewise  non-oscillatory  polynomial  recon- 
struction  of  the  point-values  of  u(x,tn)  from  their  cell  averages  {u"}  that  we  denote  by 


u(x,tn)  =  R  (x;un)  :=  X>?(aO  ■  1  ir 

3 


(2.3) 


This  reconstruction  procedure  is  at  the  heart  of  high-resolution,  non-oscillatory  central  schemes,  and  requires  the 
coefficients  of  the  polynomials  {p™(x)}  to  be  determined  so  that  R(x;un)  satisfies  the  following  three  essential 
properties. 

•  V\  —  Conservation  of  cell  averages:  p"(x)  =  it". 

•  V2  —  Accuracy:  R  (x;un)  =  u(x,tn )  +  0((Ax)r)  for  reorder  accurate  scheme,  wherever  u(-,tn)  is 
sufficiently  smooth. 

•  V3  —  Non-oscillatory  behavior  of  y]  ■  pj  ( x )  •  1  ij  which  is  characterized  in  different  ways  for  different  recon¬ 
structions.  Examples  include  TVD  reconstructions  initiated  by  van  Leer  and  Harten,  [31,  10],  Number  of 
Extrema  Diminishing  (NED)  and  shape  preserving  properties  of  Liu-Tadmor  third  order  scheme,  [22],  and 
the  higher  order  non-oscillatory  extensions  offered  by  the  Essentially  non-oscillatory  (ENO)  reconstruc¬ 
tions  of  Harten  et.  al.,  e.g.,  [9],  or  their  Weighted  versions  (WENO)  e.g.,  [27,  12],  and  their  implementation 
within  the  central  differencing  framework  e.g.,  [15,  19,  21]. 

Equipped  with  the  reconstructed  point  values  u(-,tn),  the  second  step  evaluates  the  two  time  integrals  of 
/(u(-,r).  Since  the  solution  of  (2.2)  remains  smooth  along  the  lines  {x  —  Xj}  x  [tn,tn+1],  a  simple  quadrature 
rule  suffices  to  approximate  those  integrals.  This  requires  the  intermediate  values,  {u"+,:<},  which  are  predicted 
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Xj — \  Xj ^  Xj  X j- 1_  -  1  x 

Figure  2.1:  Modified  central  differencing  of  Kurganov  and  Tadmor. 


either  by  a  Taylor  expansion  or  alternatively,  by  Runge-Kutta  solvers  of  the  ODE  uT  =  fx\x=Xj,  u(xj,  0)  =  u", 
r  >  f",  where  fx  stands  for  the  numerical  derivative  of  f{u). 

This  type  of  evolution  over  staggered  grids  based  on  cells  with  uniform  width  Ax  does  not  admit  a  semi¬ 
discrete  limit  as  At  J,  0.  Instead,  a  modified  approach  proposed  by  Kurganov  and  Tadmor  in  [17]  and  based  on 
variable  size  cells  (depicted  in  figure  2.1)  yields  the  desired  semi-discrete  limit.  Here,  we  extend  the  derivation 
of  second-  and  third-order  semi-discrete  schemes  of  [17,  15,  16],  to  schemes  of  arbitrary  order  of  accuracy  r,  [3]. 

The  first  step  passing  from  a  fully-discrete  to  a  semi-discrete  formulation,  consists  of  reducing  the  size  of  the 
staggered  cells  covering  the  Riemann  fans.  Instead  of  cells  with  uniform  width  Ax,  we  use  cells  with  variable 
width  of  order  O(At),  by  incorporating  the  maximal  local  speed  of  propagation, 


(2.4) 


here  u  ,  and  u+  ±  stand,  respectively,  for  the  reconstructed  values  of  u  to  the  left  and  right  of  the  interfacing 
3'  2  3'  2 

breakpoints,  uj+  x  :=  R(xj+ 1  ±,  un),  and  C  is  the  phase  curve  connecting  these  two  states  through  the  Riemann 
fan.  For  practical  purposes,  we  choose 

which  is  exact  for  genuinely  non-linear  and  linearly  degenerate  fields. 

This  information  allows  us  to  repartition  the  original  grid  so  as  to  distinguish  between  sub-cells  of  width 
2a",  i  At",  where  the  solution  is  dictated  by  the  non-smooth  Riemann  fan, 


~  [Xj+|,P  X0+\,r\  Kj+5  a]+  Xj+i  +  ai+\^t  ]> 

and  the  other  cells  of  width  A  Xj  :=  Ax  —  At"(a"_  i  +  a"+i ),  which  are  beyond  the  influence  of  the  Riemann 


fans  so  that  the  solution  remains  smooth, 
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In  non-smooth  regions,  the  cell  averages,  {wn_[  1 }  are  computed  via  staggered  evolution  over  the  (narrower, 

3  '  2 

O(At))  cells  covering  the  Riemann  fans, 

■•"+!  :=  fj+i'ru(x,tn+1)d. 

2  J  xn  , 


W 


3  +  i 


X  = 


3  +  l 


=  /  ^  R  (x'^n)dx-  2al  /  [/Wi"+i]/))-/WI^,i’T))](iT> 

j+^,1  •J+2 

and  in  the  smooth  regions,  direct  integration  of  (2.1)  yields  the  cell  averages  {tc"} 

:=  J-  T’  u(x,tn+1)da 
J  X71  ■, 


IX  = 


z  i  rtn+1 

f  1  R  {x\  un)dx  —  — —  /  [f(u(x™+i  (,t))  —  f(u(x™_i  r,r))]dT. 

J  xn  ,  Jtn  2  ’  2  ’ 


(2.5) 


(2.6) 


The  non- staggered  cell  averages,  },  are  then  recovered  via  a  new  non-oscillatory  piecewise  polynomial 

reconstruction  of  the  point  values  of  u(x,  tn+l)  from  the  new  cell  averages  (2.5)  and  (2.6), 


R  (x;wn+1)  =  J>”+1(x)  ■  H  +  <+1  (*) 


(2.7) 


Projection  onto  the  original  grid,  {/j},  (2.7)  yields 


«r  := 


/x  .  ,  i 

*  *i?(x; 

a? .  i 


u)n+1)dx  = 


Ax  L 


2  o7ll{x)  dx+  I  2  g”+1(x)  dx+  I  2  q™tl(x)  dx 


■Vi-*,! 


(2.8) 


Vl-*,1 


The  explicit  calculation  of  the  polynomials  <7"+1(x)  and  ^2(x)  will,  in  fact,  prove  to  be  irrelevant  for  the 
semi-discrete  formulation  that  we  seek.  As  we  shall  see  below,  only  the  cell  averages  of  these  q’s  matter,  and  in 
the  At  J,  0-limit,  these  averages  are  further  reduced  to  the  original  point  values  reconstructed  at  the  interfaces, 
p(xj±i,tn).  We  turn  to  the  details. 

Without  loss  of  generality,  we  assume  that  the  qn+1,s  are  polynomials  of  degree  r  —  1  (required  to  preserve 
the  underlying  r  -order  accuracy)  and  satisfy  properties  V\  -  V3  for  the  cell  averages  (2.5)  and  (2.6)  in  the 
corresponding  cells  Ij  and  w  We  express  these  polynomials  as  qj±l(x)  =  J2k=o^j±r(x  ~  xj± i)fcA!;  their 
coefficients  can  be  uniquely  determined  by  imposing  the  conservation  constraints1, 


/ 


^  q^(x)dx  =  w^+p,  p=  0,±i,±l,...,±^-^,  s=^’°4' 


(2.9) 


Equipped  with  this  notation,  the  first  and  last  integrals  on  the  right  of  (2.8)  read 


rx  1  r-i 


+“?_!  A<:  l  (At) 


'3  1--1 


=  (2-10) 


k= 0 


and 


Jx,  ,  1  ,  2  1 — n 


r-l  4(fe) 

r  i+i-fx-s  1  jfc+1 

t'o  L(fc  +  1)!  S+|) 


r— 1 


Jx.  ,  1  — an  ,  At  - 

j+i  3+%  k—0 


=  E  <2 


(fc  +  !)!  K  3+h’ 


1The  choice  of  symmetric  intervals  around  Ij  yield  methods  of  odd  order.  For  methods  of  even  order,  it  suffices  to  sacrifice  ■ 


of  the  side  constraints,  p  =  f 
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For  the  second  integral  in  (2.8),  we  observe  that  by  property  V\  the  polynomial  qrt+l{x)  must  satisfy 


++1( x)  dx  =  A XjWj 


n+1 


=  [Ax  -  At(a]_l  +  a]+1_)\w] 


-n+1 


(2.12) 


Substituting  (2.10)-(2.12)  into  (2.8),  we  arrive  at  a  fully-discrete,  non-staggered  approximation  (of  arbitrary 
order  of  accuracy  r)  of  the  cell  averages  in  (2.2), 


1 - 1 


-n+1 


=  £ 

k—0 


A(A  t) 


k  r 


(k  +  1y. iK-d‘+I++(-1)',K+d‘+1+]  +  [1-A(»i-5+“«i>]<1'  A:=  +  <2-13> 


Given  the  appropriate  reconstructions  R(x;un ),  R(x',uin),  and  quadrature  rules  to  evaluate  the  flux  integrals 
in  (2.5)  and  (2.6),  then  (2.13)  constitutes  a  family  of  fully-discrete  non-oscillatory  central  schemes,  e.g.,  [11, 
17].  This  type  of  schemes  admit  a  considerably  simpler  and  more  versatile  semi-discrete  limit,  dUj/dt  := 
liniAt^o(%(t"  +  A t)  —  Uj(tn))/At.  To  this  end,  we  rewrite  (2.13)  as 


duj 

dt 


=  lim 


1 


J  +  +  T&«r  -  «<)  -  +z(»U  +  «"*i) 


Ai— >0  1  Aa’  3  2  i- 


,(0) 

"3+^  3+ 


1 


-n+1 


At 


r—  1 

+£ 

k=l 


A(A  t) 


k—  1  r 


(Jfe  +  1)!  L  3~i 


(ati)k+lA?-l  +  (-l)fc(++i)fc+1£+i 


3  +  i 


=  lim 


1 


+  —  «+i  -  <)  -  — (a+,  +  a"+l)w 


l(0) 


1 


At^O  1  Ax K  3-2  3-2  3+2  3+2 


At 


,n+l 
3 


Ax  3  2 


1 


3+1 


—  n+1 


A:r  3  2 


o+i 


(2.14) 


Since  the  width  of  the  local  Riemann  fans  approaches  zero  in  the  At  |  0-limit,  i.e.,  a;”+i  (,i"+i  r 
obtain  by  conservation  of  cell  averages  property  V\,  A ( L  =  wnfl  so  that  (2.14)  yields, 

3± 2  3± 2 


dun 

=  lim 


dt 


*=.  i  h<aU*?l +  +  ItW  - »?)  -  h{aU  + 


}■ 


X3+b  We 


(2.15) 


As  At  |  0  we  end  up  with  limiting  values  which  depend  solely  on  the  reconstructed  values  R(x,  u )  on  both 
sides  of  the  interfacing  breakpoints,  x.J±  i ,  given  by 


+++!,+)  -  Pj(xj+i)  =:  «j+i(<) 

“(*"+*,»■>*)  -  p"+i  +  i)=:«+(t); 


(2.16) 


here  p”(a;)  and  p"+1(a:)  are  the  polynomial  reconstructions  of  R{x\un )  at  the  original  time  step  tn  introduced 
in  (2.3).  The  corresponding  averages  in  (2.5)  approach  wn+  ~  (A,/++ 1  /2a”+1  where  (A/)J+i  :=  /(w++i)  — 


/(u  i),  while  the  smoother  averages  in  (2.6)  are  of  order  t++1  ~  tt"  —  Af(A/)J+i .  Inserted  into  (2.15),  we 

end  up  with  a  semi-discrete  limit  which  can  be  written  in  its  final  conservative  form, 


duj 

dt 

here  -fi+i  (f)  are  the  numerical  fluxes  given  by 


ff#+4(i)-JWt) 


Aa’ 


(2.17) 


+++*)  “ 


f(ut+±  (*))  +  /(++!  W)  a,-+i  (t) 


i+2 


+7+3 


-++!+]• 


J+5 


(2.18) 


Once  again,  we  emphasize  that  among  its  several  advantages,  the  versatility  of  this  semi-discrete  formulation 
is  the  independence  of  the  specific  reconstruction,  (2.3),  and  evolution  algorithms  which  are  utilized  to  recover 
the  interface  values,  u~+ 1  and  and  to  evolve  the  cell  averages  in  (2.17).  A  description  of  our  particular 
choices  for  these  algorithms  is  provided  in  §3. 
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Figure  2.2:  Floor  pan  for  modified  central  differencing  in  two-space  dimensions. 


2.2  Two-dimensional  Schemes 

Starting  with  a  general  hyperbolic  conservation  law  in  two  space  dimensions, 

ut  +  f{u)x  +g{u)z  =  0,  (2.19) 

we  proceed  as  in  §2.1,  and  consider  the  sliding  averages  of  u  over  the  cells  Ijtk  :=  [xj- 1/2,  £.7+1/2]  x  [zk_ 1/2,  Zfc+1/2], 


u™  k  :=  T  u(x,z,tn)  dx  dz, 

J  1, 1. 


where  Ax  and  A;:  represent  the  space  scales  in  the  x-  and  ^-direction  respectively. 

As  in  the  one-dimensional  case,  we  seek  a  semi-discrete  formulation  for  the  evolution  of  the  cell  averages 
of  the  conservation  law  (2.19).  We  begin  by  incorporating  the  information  provided  by  the  local  speeds  of 
propagation,  that  we  approximate  by 

a”+iife:=max{p(|^(«^life)),p(|^(^fc))},  :=  max  {d(ff(<4+i))  ,  p(^«fc))  };  (2-20) 

where  the  cell  interface  values  in  the  z-  and  x-direction, 


ufk  '■=  Pj,k(xj’zk+±),  uj,k  '■=  Pj,k(xj>zk- i). 


',k  '  Pj,kiXj+  iizk')j 


J-Yk  ~Pj,k(xj-^zk) 


(2.21) 


are  calculated  via  a  non-oscillatory  piecewise  polynomial  reconstruction, 

R{x,z-,un)  =  J2plk(x>z)  (2-22) 

j,k 

the  polynomials  {prJk(x,  z)}  are  determined  so  that  R(x,  z; «")  satisfies  properties  analogous  to  V\  -  V3  above. 
This  information  allows  us  to  separate  between  regions  of  smoothness,  depicted  as  D]k  in  Fig.  2.2,  and  regions 
of  non-smoothness  ,  depicted  as  the  shaded  regions.  Nine  sets  of  cell  averages  are  calculated  as  follows: 

•  In  the  clear  shaded  regions,  DJ±  1  the  solution  is  not  smooth  in  both  directions;  staggered  evolution 
over  the  control  volumes  Dj±i  k±i  x  [tn,tn+1]  is  used  to  obtain  the  new  cell  averages  {^"^"1  fc±i}, 

•  In  the  dark  shaded  regions,  Dj±  i  k  and  Dj  k±  1,  the  solution  is  smooth  only  in  one  direction;  staggered 
evolution  is  used  along  the  non-smooth  interfaces  to  obtain,  respectively,  the  cell  averages  {to?_£i  k}  and 

•  To  evolve  the  smooth  part  of  the  solution,  we  integrate  the  polynomial  p™k(x,z)  and  the  corresponding 
fluxes  over  the  non  rectangular  control  volume  Dj  k  x  [tn,  t"+1];  these  cell  averages  are  denoted  by 
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This  approach  allows  us  to  form  a  new  polynomial,  denoted  R(x,  z\  utn+1),  which  is  reconstructed  from  these 
smooth  and  non-smooth  portions  of  the  solution,  and  to  reproject  it  back  onto  the  original  grid-cells, 


= 


J-  R  (x,  z;wn+1)  dx  dz. 
J  I.,.* 


(2.23) 


The  resulting  non-staggered  fully-discrete  scheme  admits  a  more  versatile  semi-discrete  limit  as  At  — >  0, 

d_...  - Hu, *(*) 


Ax 


A^ 


(2.24) 


with  numerical  fluxes 
1 

12  L' 


= 


+  / (uf,k  W)  +4 if{Uj  +  i'k(t))  +  +  /(wf+i,fcW)  + 

(Uj+T,k  W  -  Uj'fit)  +  Muf+Uk (t)  -  ufik{t))  +  «?£,*(*)  -  u? f  (f))  ,  (2.25) 


S+|,fcW 
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3,k+i 


(t)  =  [gKMt))  +  <?KT  W)  +  ^Kk+1  (t))  +  <?(«&(*)))  +  <7(^|+1  (*))  +  s(<n*)) 


bj,*+i  (*) 


12- —  (wiA+i(*)  “  (*)  +  4(uffc+iW  -  +  uj,k+i(t)  -  (*))  •  (2-26) 

This  particular  version  of  the  numerical  fluxes  results  from  using  Simpson’s  quadrature  rule  to  approximate 
the  integrals  of  the  fluxes  /  and  g  along  the  cell  boundaries  [zk_i,zk+ 1]  and  [Xj_i,Xj+i]  respectively,  and  it 
incorporates  information  from  the  corner  interface  values, 


,NE 


,k  •  )?  ^j,k  ’  i  %k—  ^ )  > 

•Sk  ~Plk(xj+^,zk_i),  ufj?  :=  Pjtk(Xj_i,  zk+i), 


(2.27) 


into  the  scheme.  These  corner  values  are  recovered  via  the  non-oscillatory  reconstruction  R(x,  z ;  un)  =  XI  Pj,k  (x;  z)- 
1  /  fe,  which  may  coincide  with  the  original  reconstruction,  R{x ,  ft11),  or,  alternatively,  interpolate  the  cell  val¬ 

ues  along  the  diagonal  directions  shown  in  figure  2.3(b)  so  as  to  prevent  the  onset  of  spurious  oscillations  along 
those  axis. 


(a) 


Figure  2.3:  (a)  reconstruction  in  x-  and  z-  directions.  (6)  diagonal  axis  for  reconstruction. 


3  Implementation  of  Semi-discrete  Central  Schemes 

The  framework  of  semi-discrete  scheme  described  above  entertains  a  wide  range  of  options  for  the  implementation 
of  their  two  main  ingredients:  non-oscillatory  reconstruction  and  evolution.  In  this  section  we  provide  some 
examples  of  non-oscillatory  reconstructions,  evolution  routines  and  magnetic  field  corrector  that  we  implemented 
for  computing  the  solutions  of  the  system  (1.1)  -  (1.4)  presented  in  §4  and  §5  below. 


3.1  Third-order  CWENO  Reconstruction 


Our  first  choice  for  the  reconstruction  of  the  point  values  of  u  is  the  third-order  CWENO  polynomial  recon¬ 
struction  of  Kurganov  and  Levy,  [15].  A  piecewise  quadratic  polynomial  that  satisfies  the  essential  properties 
Pi,  P2,  and  P3  above  is  formed  as  follows:  in  each  cell  Ij  =  [x-_i,a;-+i],  the  polynomials  {pj( x)}  in  (2.3)  are 
written  as  a  convex  combination  of  three  polynomials  PL(x),  Pc(x),  and  PR(x), 

Pj{x)  =  wLPR(x)  +  wcPc(  x)  +  wRPR(x),  ^2  wi  =  l,  (3.1) 

ig{L,c,R} 


where  the  linear  polynomials 


U ^ 

PL(x)  =■»’•'  +  3  J~1(x-Xj))  and  PR(x)  =  u™  +  - 3-(x-Xj), 


Ax  v  J  Ax 

conserve  the  pair  of  cell  averages  u”_1,  u”  and  w",  u”+1  respectively,  and  the  parabola  centered  around  Xj , 


(3.2) 


Pc(x)  =  Uj  —  -^(«"_i  —  2m"  +  Uj+i) 


"i+1  “^-(x  -  Xj)  +  UU  || +  +1  (x  -  Xj)2, 


2Ax 


is  determined  so  as  to  satisfy 


(3.3) 


Cl-Pl(x)  +  CrPr(x)  +  (1  -  Cl  -  Cr)Pc(x)  =  Uj  +  u'j  (x  -  Xj)  +  ~u”(x  -  Xj)2, 
where  u" ,  it'  and  u"  approximate  the  point  values  u(xj,tn),  ux(xj,tn ),  and  uxx(xj,tn)  given  respectively  by 


””  :=  a?  -2a?  +  «“«>• 

7/n  _  j-.n  -n  _  Qy.n  .  y.n 

«3+i  and  »  1 

3  2Ax  3  Ax2 


(3.4) 


We  guarantee  the  conservation  of  the  cell  averages  Uj_1,  u"  and  u"+1  by  any  symmetric  choice  of  the  weights 
Cj  (e.g.,  cL  =  cR  =  1/4,  cc  =  1/2),  and  the  accuracy  property,  P2,  in  smooth  regions,  [20] 

Pj(x)  =  u(x,tn)  +  0((Ax)3),  xG  [ar^i.a^+i]. 

The  non-oscillatory  behavior  of  { p"(x)}  property  P3,  is  attained  with  the  non-linear  weights 


Em  am 


with  a*  = 


(e  + 1  Si)* 


i,  m  G  {l,  c,  r}, 


(3.5) 


where  t«l  prevents  the  denominator  from  vanishing  (for  the  calculations  in  §4  and  §5  we  choose  e  =  10  6), 
and  the  smoothness  indicators,  IS),  provide  a  local  measure  of  the  derivatives  of  Pj(x), 


I  Si  =  ^2  f  (Ax)21  1(P?\x))2dx,  i  G  {l,  c,  r},  (3.6) 

/= 1  ’’ 1  3 

switching  automatically  to  the  second-order  reconstructions  PL  and  PR  in  the  presence  of  steep  gradients  and 
avoiding  the  onset  of  spurious  oscillations,  [12,  15,  20].  In  this  case,  they  read 

ISL  =  (u]  -  u^)2,  ISR  =  («”+1  -  u”)2, 

I  So  =  y  («?+ 1  -  2  u]  +  u -_,)2  +  \(u]+1  -  u $U)2-  (3.7) 


Remarks: 
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1.  For  this  particular  reconstruction,  the  one-dimensional  interface  values,  u  .  j  (t)  and  u+  x  (<),  required  in 

3'  2  3\  2 

the  numerical  flux  (2.18)  take  the  explicit  form 


1 


j+ 


I  (tn)  Pj(Xj+ 1)  —  -wL(3u"  —  Uj_  i)  +  —  wc(— 5m"+1  +  8  uj  —  u”_r)  +  -rcR(u"+1  +  «”), 


12 

1 


J+2 


i  (*")  :=  p”+i(^+i)  =  ^wL(w"+1  +  «")  +  — Wc(-w"+2  +  8u"+1  -  5Uj)  +  -%(-mJ+2  +  3u"+1). 


2.  In  the  case  of  systems  of  equations,  the  choice  of  these  smoothness  indicators  is,  indeed,  non-trivial. 
The  different  conserved  quantities  involved  in  any  particular  system  (density,  momentum,  etc.)  may 
very  well  develop  discontinuities  at  different  points  throughout  the  solution  domain  that  do  not  affect 
every  other  quantity  in  the  system.  Hence,  choosing  the  smoothness  indicators  individually  for  each 
quantity  may  cause  the  scheme  to  use  different  stencils  for  different  quantities  during  the  same  evolution 
step.  This  can  be  avoided  by  either  using  global  smoothness  indicators,  e.g.,  an  (norm-scaled)  average  of 
the  individual  ones,  or  by  identifying  those  that  are  physically  relevant,  such  as,  e.g.,  the  density  near  a 
contact  discontinuity  (see  §4  for  our  particular  choices  in  the  one-dimensional  case  and  [19]  for  an  in-depth 
discussion  about  the  selection  of  smoothness  indicators  for  systems  of  equations). 

In  two-space  dimensions,  one  possibility  for  calculating  the  interface  values  in  (2.21)  and  (2.27)  is  to  apply 
this  one-dimensional  third-order  reconstruction  dimension-by- dimension  in  the  x-  and  z-  directions  (Fig.  2.3 
(a)),  and,  if  so  desired,  in  the  two  diagonal  directions  of  the  coordinate  frame  displayed  in  figure  2.3(b).  In 
two  dimensions,  however,  we  note  that  the  constant  term  in  the  central  parabola,  (3.3),  requires  an  additional 
correction  (in  the  transverse  direction  of  the  reconstruction)  in  order  to  guarantee  third-order  accuracy  when 
the  point  value  u(xj,  Zk,tn)  is  recovered  from  the  neighboring  cell  averages, 


u(xj,zk,  tn)  =  M  -  2 u\k  +  u™+  i>fc)  -  —  («]•>_ i  -  2 ulk  +  u£fc+1)  +  e>((max{Aa’,  A^})^. 
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That  is,  along  the  line  z  =  zk,  the  reconstruction  in  the  ^-direction  is  carried  out  using  the  linear  polynomials 
PL(a :,  Zk)  and  PR(x,  zk)  as  in  (3.2),  and  the  parabola  Pc{x,  zk),  given  by 


Pc(x,Zk)  —  llj'k  —  —  (Uj_lk  —  2u™k  +  U^+lk)  —  —  («”,*_!  —  2w™fc  +  u™k+i) 


“j+l.fc 


k~Uj- l,k,  s  -  2«"  +  Uj+l,k  ,  -,2 

2Aa‘ 3  (X  ~  X’]  +  - Ax2  (X  (3’8) 

e  same  reconstruction  is  carried  out  in  the  ^-direction,  holding  x  =  Xj  fixed. 

The  corresponding  polynomials  in  the  diagonal  directions  are  given  by  (consult  [16]  for  additional  details) 


with 


PL{x,z) 

Pc{x,z) 


p™k{x,  z)  =  wLPL{x,  z)  +  wcPc{x,  z)  +  wRPR(x,  z), 

•** + - - Zk)  +  2a5(i  -  *>>) 

1 


(3.9) 


ui,fc  ~  J^(uj'+i,k+i  —  +  —  12  ~  2u")fe  +  w"+i,fe_i) 

,  “t+i.fc+i  _  /A  A 

+ - 5a -  5a5<"  -  +  5 X5(x  ~ 


2A  V2Azv  ^  '  2Axk 


Pr(x, z) 


A2 

Uj+ i,fc+i  _  Uj,k  (  A 


A 


A 

VtAz Ci:2Ar 


{x  - 


h 


(3.10) 


10 


along  the  SW-NE  axis,  and 


PL(x,z) 

Pc{x,z ) 


4  Or,  z) 


xj,k  ' 


Uj,k  Uj+l,k-l 


~ Zh)  ~  ~ Xi)) 


A 


A 


Uj,fc  ~  Y2^U^+1’k+1  ~  ^UIk  +  U]-l,k-l)  ~  J^(ulj-l,k+l  —  2ll™k  +  u^+lk-i) 


Uj-l,k+l  Uj  +  l,k-l 

2A 

T:n  —  9fin  4-  fin 

“j-l.fc+l  LLj,k  t  »-,  +  !, fc-1 


A2 


J/j,k  ' 


Mi+l,fe+l 


-  <k  (  A 


A 


A 


(3.11) 


along  the  SE-NW  axis  where  A  =  (Acc)2  +  (A z)2. 

The  non-linear  weights,  u1,,  are  also  calculated  direction  by  direction  according  to  the  one-dimensional  recipe 
(3.5)  with  the  smoothness  indicators  given  by  the  norm-scaled  average  of  the  componentwise  indicators,  IS\m\ 
calculated  as  in  (3.7), 


1 


d 

/s-  =  sE 

m=  1 


,(m) 


2  +  e 


IS 


(m) 


i  e  {l,  c,  r} , 


(3.12) 


where  u ^  stands  for  the  mth  component  of  u,  and  ||ibm)||2  =  J2j  k  \upk  |2  AxAz  represents  its  1 2  norm  over 
the  solution  domain. 


3.2  CWENO  Fourth-order  Reconstruction 

In  §5  below  we  demonstrate  the  versatility  of  the  central  semi-discrete  scheme  (2.24)-(2.26)  by  implementing 
the  dimension-by-dimension  third-order  reconstruction  described  above  and  the  genuinely  multi- dimensional 
fourth-order  CWENO  reconstruction  of  Levy  et.  al.,  [21,  3],  that  recovers  the  interface  point  values  in  (2.21) 
and  (2.27)  via  a  bi-quadratic  polynomial  interpolant  satisfying  the  desired  conservation,  accuracy  and  non- 
oscillatory  properties  in  two-space  dimensions.  Following  is  an  outline  of  this  fourth-order  reconstruction, 
originally  developed  by  Levy  et.  al.  within  the  central,  fully-discrete,  staggered  framework. 

In  each  cell  Ij,k,  the  piecewise  polynomial  interpolants  in  (2.22),  p™k(x,z),  are  written  as  a  convex  combi¬ 
nation  of  nine  bi-quadratic  polynomials,  Pj+p^+q{x,  z),  that  conserve  the  cell  averages  {u"+  fe+  and 

approximate  the  pointvalues  of  u(x,z,tn )  within  fourth-order  accuracy, 


Pj,k(x,z ) 


1 

W^’kPj+P,k+q(x,  z), 

p,q--l 


1 


p<q=- 1 


W 


P,< 7 

j,k 


>0; 


(3.13) 


where  the  polynomials  Pjpjx,  z)  have  the  form  (omitting,  for  simplicity,  the  j  and  k  indeces  in  the  coefficients 

{bm}), 


Pj,k{x,z)  =  b0+  bi{x-Xj )  +  b2(z-zk)+  b3(x  -  Xj)(z  -  zk) 
+  64(21  -  xj)2  +  b5(z  -  zk)2  +  b6(x  -  Xj)2 (z  -  zk) 
+  b7(x  -  Xj){z  -  zk)2  +  bs(x  -  Xj)2(z  -  Zk)2- 

The  9  conservation  constraints, 

[X.+  i+pAx  rz^i+qAz 

f  f  PjAx,  z)dzdx  =  Uj+p  k+q,  p,q=- 1,0,1, 

J  a:._  1 +pA xJ  zk_  i_ -\-qAz 


(3.14) 


(3.15) 
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uniquely  determine  the  9  coefficients  {6m}, 


Ax2 


Az2 


- Urr - UZz 

24  xx  24 


Az2  „ 

Cl  =  Ux  “  ~  UXzzi 


h  = 

^5  = 

b7  = 

where  the  divided  differences 


24 


1  „  Aa’2  „ 

~Z^ZZ  .  Q  XXZZl 

Z  4o 


Wj+l,/c  u1j-l,k 


Ax2  A z2  , 

24^  U'XXZZ’) 

Ax2  „ 

^2  =  ^2:  24  'Ujxxz’> 

l ,  A22 „ 

^4  —  "Z^xx  A  0  'Ujxxzzi 

Z  4o 

^6  —  ’^V'xxzi 
^8  —  ^'^'31212:2:5 


^fc+l 


*xxj!k 


J'xzj,k 


t'XXZj'k 


*XZZjjk 


2Ax 

r,n  _  9 7-,n  .  z-.n 

xj+l,k  Laj,k  ^  uj-l,k 


2Az 


ulk+ 1  -  2 «",fc  + 


Ax2 


Az2 


'//n 

Uj+l,k+l 


ij+ i,fc— 1  Mj-i,fc+i 


4A:rAz 

(“j+i.fe+i  —  2-^i,fc+i  +  ^j-i,fc+i)  —  (^yfi,fc-i  —  2^i,fc- 1  + 

2Aa;2Az  ’ 

(“i+i.fe+i  —  2^yfi  ,fc  +  ^yfi.fc-i)  —  (^i-i,fe+i  —  2^i-i,fc  +  ^i-i,fe-i) 


1 


^XXZZj^k 


2AxA  z2 

Ax2A~2  [^i+1>fe+1  ~~  2^a+i,fc  +  ^i+i,fc-i)  ~  2C“j,fc+i  —  2“i,fc  +  Uj,k- 1) 

+(^j-i,fc+i  ~  2“i-i ,fc  +  ^j-i,fe-,i)]  1 

serve  to  approximate  the  values  of  u  and  its  partial  derivatives  at  the  points  (a ’j,Zk)  within  the  fourth-order 
accuracy  constraints  of  the  method. 

As  in  the  third-order  reconstruction,  the  non-linear  weights  in  (3.13),  wp’^,  are  computed  so  as  to  provide 
maximum  accuracy  in  smooth  regions  and  prevent  oscillations  in  the  non-smooth  regions  by  eliminating  the 
contribution  of  polynomials  with  steep  gradients  across  the  cell  interfaces.  For  each  cell  nine  weights  are 
calculated, 


wpi  = 

J,k 


ryP’i 

aj,k 


with  a?’i  = 


„P>4 

cj,k 


m,n,p,q  =  -1,0,1, 


.  (e+is$)»- 

and  the  linear  coefficients,  c?’^,  chosen  so  that  symmetry  guarantees  fourth-order  accuracy, 


(3.16) 


„o,o  ^  _ 

cj,k  ~  2 


and 


=  ->  Vp,9^0; 


16 


convexity,  ^  4’fe  =  1>  cj’I  —  guarantees  conservation. 


Global  smoothness  indicators  are  chosen  according  to  (3.12).  Unlike  the  third-order  reconstruction  where 
the  exact  calculation  of  the  individual  smoothness  indicators  is  trivial  and  renders  simple  formulas,  i.e.,  (3.7), 
in  this  fourth-order  case,  the  exact  calculation  of  the  local  L2-norm  of  the  partial  derivatives  of  the  polynomials 
Pp{!  is  impractical.  Instead,  a  Gaussian  quadrature  formula  is  used  to  approximate  the  integrals 


IS? fc  =  /  (\dxPj+Ptk+q\2  +  \dzPj+p,k+q\2  +  {Ax)2\dxxPj+Pik+q\2  +  ( Az)2\dzzPj+Ptk+q\ 2)  dx  dz. 

Jij.k 

This  approach  still  serves  the  purpose  of  automatically  detecting  and  re-directing  the  numerical  stencils  in  the 
direction  of  smoothness,  thus  preventing  the  onset  of  oscillations. 
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3.3  Time  Evolution:  SSP  Runge-Kutta  Solvers 


With  the  interface  values  (2.16),  (2.21),  and  (2.27)  computed  with  either  one  of  the  above  CWENO  recon¬ 
structions,  we  denote,  for  any  grid  function  w  =  {uij(t)},  the  numerical  fluxes  on  the  right  of  (2.17)  and  (2.24) 
by 


C[w{t)] 


Hj+i(w(t))  -Hj_i{w{t)) 

Ax 


(3.17) 


and 


C[w(t)}  =  - 


HJ+i 


,k 


(w(t))  -  H*_; 


,k 


( w(t )) 


if* 

j,k+- 


-  Hzj  k  l{w{t)) 


Ax 


A  z 


(3.18) 


respectively,  and  evolve  the  solution  u(0)  :=  u"  from  tn,  to  tn+1  with  an  appropriate  ODE  solver.  For  the 
numerical  calculations  reported  in  §4  and  §5,  we  choose  the  Strong  Stability-preserving  (SSP)  Runge-Kutta 
discretizations  [26,  8].  In  particular,  for  the  third-order  results,  we  evolve  the  solution  according  to 


uW  =  u(o)  +AtC[M(0)], 

m(2)  =  U(1)  +  ^(_3C[U(0>]  + 

Un+1  u(3)  =  u(2)  +  ^(_(7[u(0)]  _  cr[u(i)]  +  8  c[u(2)}).  (3.19) 

and  for  the  fourth-order  results  of  §5.2 

u*1)  =  u<°>  +  ^-C[uW], 

uS 2)  =  «(1)  +  y(-c[«(0)]  +  c[m(1)]), 
u(3)  =  u(2)  +  ^*(_cf[u(i)]  +  2  C[u(2)]), 

un+ 1  :=  u(4)=w(3)  +  — (C[u(0)]+2C[u(1)]-4C[m(2)]+C[u(3)])  (3.20) 
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is  used  to  evolved  the  values  reconstructed  via  the  fourth-order  CWENO  interpolation  algorithm  of  Levy  et. 
al.  of  §3.2. 


3.4  Magnetic  Field  corrector 

In  addition  to  reconstruction  and  evolution,  our  central  schemes  requires  the  correction  of  the  magnetic  field, 
B,  in  order  to  guarantee  the  solenoidal  constraint  V  •  B  =  0.  We  satisfy  this  constraint  by  updating  the  cell 
averages  of  the  magnetic  field  at  the  end  of  each  time  step,  replacing  the  computed  B,  with  its  divergence-free 
projection,  Bc.  To  this  end  the  so  called  Leray  projection  is  carried  out  by  solving  the  Poisson  equation 

A0=-V-B  (3.21) 

with  the  appropriate  boundary  conditions.  Since  the  coplanar  structure  of  the  problem  guarantees  dBy/dy  =  0, 
only  the  components  Bx  and  Bz  need  to  be  updated.  We  use  a  fast  Poisson  solver  for  the  standard  five  point 
discretization  of  the  potential  (f>,  and  central  differences  for  the  divergence  of  the  magnetic  field,  V  •  B.  The 
corrected  -  divergence  free  —  magnetic  field,  Bc,  is  then  realized  as 

BC  =  B  +  V<^>.  (3.22) 

Applying  the  divergence  operator  V  =  (d/dx,  d/dz)  to  both  sides  of  (3.22),  one  can  easily  verify  that  V-Bc  =  0. 
We  use  here  only  one  out  of  several  methods  to  enforce  the  so-called  ‘constrained  transport’  and  we  refer  the 
reader  to  [5] ,  [25]  and  the  references  therein  for  a  general  discussion  and  [30]  for  handling  the  solenoidal  constraint 
in  the  context  of  MHD  schemes  over  staggered  grids. 

The  solution  of  (3.21)  is  computed  with  a  C  translation  of  Fishpack,  a  collection  of  Fortran  subprograms 
which  utilize  cyclic  reduction  to  directly  solve  second-  and  fourth-order  finite  difference  approximations  to 
separable  elliptic  Partial  Differential  Equations,  [1,  2]. 
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4  One-dimensional  Numerical  Results 


In  one  space  dimension,  equations  (1.1)— (1.4)  admit  the  conservative  form  (2.1)  with 


u  =  (p,pvx,pvy,pvz,By,Bz,e)T , 

f(u)  =  (pvx,  pvl  +  p*  -  Bl,pvxvy  -  BxBy,  pvxvz  -  BXBZ,  Byvx  -  Bxvy, 
Bzvx  BxvZi  (c  ~t -  p  'jVx  Bx{Bxvx  T  Byvy  Bzvz')')  , 


(4.1) 

(4.2) 


where  p*  =  p+  ^B2  stands  for  the  total  pressure  (static  plus  magnetic). 

In  this  section  we  present  numerical  simulations  of  this  one-dimensional  MHD  equations,  (2.1)-(4.2).  The 
results  were  obtained  using  the  semi-discrete  central  scheme  (2.17)-(2.18)  with  the  pointvalues  u^+l(t)  and 


u  i  (f),  (2.16),  computed  via  the  CWENO  third-order  reconstruction  of  §3.1,  and  evolved  according  to  Shu’s 
.2 

third-order  SSP  Runge-Kutta  scheme  as  outlined  in  §3.3.  The  schemes  are  implemented  for  computing  the 
approximate  solution  of  two  coplanar  shock  tube  MHD  models  described  by  Brio  and  Wu  in  [6].  We  use  a 
uniform  grid  in  the  space  discretization,  and  in  both  cases  we  choose  the  time  step  dynamically  with  CFL 
restriction 


At  =  0.475 


Ax 

maxfc  |  a,k  (u)  \ ' 


(4.3) 


where  {07- (it)}  a-  are  the  eigenvalues  of  the  Jacobian  matrix  of  /(it). 


4.1  Brio— Wu  Shock  Tube  Problem 


The  first  one  dimensional  Riemann  problem  we  consider  consists  of  a  shock  tube  with  two  initial  equilibrium 
states,  ui  and  ur,  given  by 


I  (1.0, 0,0,0, 1.0, 0,1.0)T  for  a;  <  0, 
\  (0.125,  0,0,0, -1.0,  0,0. 1)T  for  a:  >  0, 


(4.4) 


and  complemented  with  the  constant  values  of  Bx  =  0.75  and  7  =  2.  The  problem  is  solved  for  x  £  [—1, 1] 
with  800  grid  points,  and  numerical  results  are  presented  at  t  =  0.2.  Figure  4.1  show  the  density,  the  x —  and 
y— velocity  components,  the  y— magnetic  field,  and  pressure  profiles. 

The  hydrodynamical  data  of  this  problem  is  the  same  as  that  in  Sod’s  shock  tube  problem  of  gas  dynamics. 
The  variety  of  MHD  waves,  however,  poses  a  considerable  challenge  for  high-resolution  such  as  the  black-box 
central  schemes  described  in  this  paper.  The  solution  of  this  problem  consists  of  a  left-moving  fast  rarefaction 
wave  (FR),  a  slow  compound  wave  (SM)  which  results  from  an  intermediate  shock  that  changes  By  from  0.58 
to  —0.31  and  a  slow  rarefaction  that  changes  By  from  —0.31  —0.53,  a  contact  discontinuity  (C),  a  right-moving 
slow  shock  (SS),  and  a  right-moving  fast  rarefaction  wave  (FR).  Note  that  the  solution  to  this  problem  is  not 
unique  if  Bz  and  vz  are  not  identically  zero. 

The  results  in  figure  4.1  are  in  agreement  with  those  previously  reported  in  [4],  and  are  comparable  with 
the  second  order  upwind  computations  of  Brio  and  Wu  in  [6],  and  with  the  fifth  order  WENO  computations 
presented  by  Jiang  and  Wu  in  [14].  In  fact,  the  present  results  show  a  better  control  of  the  oscillations  that 
typically  appear  at  the  trailing  edge  of  the  right-moving  fast  rarefaction  wave  when  high  resolution  schemes 
are  employed.  This  better  better  control  is  due  to  our  choice  of  smoothness  indicators.  Our  various  numerical 
experiments  show  that  -in  the  absence  of  characteristic  information-  the  average  of  the  smoothness  indicators  of 
the  density,  p,  and  transverse  magnetic  field,  By,  (scaled  by  their  I2  norm)  is  the  best  combination  for  selecting 
a  single  stencil  for  reconstruction  and  evolution  of  all  the  conserved  quantities  in  the  system.  Our  numerical 
confirm  the  ability  of  central  schemes,  whether  in  their  fully-discrete  or  semi-discrete  formulations,  to  capture 
the  main  features  of  the  discontinuous  MHD  solutions,  while  avoiding  any  characteristic  information  other  than 
an  estimate  of  the  maximal  speed  of  propagation,  max*,  |afc(w)|. 


4.2  Brio  Wu  High  Mach  Shock  Tube  Problem 

The  following  shock  tube  model  proposed  by  Brio  and  Wu  in  [6] ,  is  commonly  used  to  check  the  robustness  of 
the  numerical  schemes  for  high  Mach  number  problems.  The  initial  equilibrium  states,  ui  and  ur ,  are  given  by 


I  (1.0, 0,0,0, 1.0,  0,1000)T  for  2;  <  0, 
\  (0.125,  0,0,0, -1.0,  0,0. 1)T  for  2:  >  0, 


(4.5) 
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Figure  4.1:  Results  of  Brio-Wu  shock  tube  problem  at  t  =  0.2  computed  with  800  grid  points  using  third-order 
CWENO  reconstruction,  (3.1),  and  Shu’s  SSP  Runge-Kutta  solver,  (3.19). 


and  complemented  with  the  values  of  Bx  =  0,  and  7  =  2.  The  Mach  number  of  the  right-moving  shock  wave  is 
15.5.  If  the  plasma  pressure  is  replaced  by  the  sum  of  the  static  and  magnetic  pressures-denoted  by  p*  above, 
this  model  becomes  a  standard  hydrodynamical  Riemann  problem.  The  solution  is  presented  at  t  =  .012, 
x  £  [—1, 1],  with  200  grid  points  and  with  CFL  number  0.475,  consult  (4.3). 

The  solution  of  this  second  Riemann  problem  consists  of  a  left-moving  fast  rarefaction  wave  (FR),  followed 
by  a  tangential  discontinuity  (TD),  and  a  right  moving  fast  shock  (FS)  with  Mach  number  15.5.  Across  the 
tangential  discontinuity,  the  density,  the  magnetic  field  and  the  pressure  can  change,  but  both,  the  fluid  velocity 
and  the  total  pressure,  p+  4^  are  continuous. 

As  in  the  previous  problem,  our  results  are  comparable  to  those  reported  in  [4] ,  [6] ,  and  [14]  with  second-order 
central  and  second-  and  fifth-order  upwind  scheme,  and  demonstrate  the  robustness  of  the  schemes  described 
in  §2. 


5  Two  Dimensional  Numerical  Results 

In  two  space  dimensions,  equations  (1.1)-(1.4)  admit  the  conservative  form  (2.19)  with 

U  (^P,pVx,pVy,pVZ,Bx,I3y,I3z,C)  ,  (5.1) 

/(u)  ( pry ,  pvx  -I -  p  Bx ,  pvxVy  BxBy ,  pvxvz  Bx B z ,  0,  ByVx  BxVy , 

Bzvx  Bxvz,  (e  +  p  )ux  Bx(Bxvx  -(-  ByVy  -t-  Bzvz) )  ,  (5.2) 

g(u)  =  {pvz,pvzvx  -  BzBx,pvzvy  -  BzBy,pv2z  +p*  -  B\,Bxvz  -  Bzvx, 

Bvvz  -  Bzvy,  0,  (e  +  p*)vz  -  Bz(Bxvx  +  Byvy  +  Bzvz))T .  (5.3) 

In  this  section  we  present  the  numerical  solutions  of  two  prototype  models  of  two-dimensional  MHD  equa¬ 
tions.  For  the  first  problem  — the  Kelvin-Helmholtz  instability  with  transverse  magnetic  field  configuration, 
we  consider  two  different  sets  of  boundary  conditions  in  the  x— direction:  periodic  in  the  first  case  and  a  free 
outflow  boundary  in  the  second  convective  setup.  The  second  problem  was  introduced  by  Orszag  and  Tang 
in  [24]  as  a  simple  model  to  study  MHD  turbulence,  and  has  become  a  standard  model  to  validate  numerical 
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Figure  4.2:  Results  of  Brio- Wu  high  Mach  problem  at  t  =  0.012  computed  with  200  grid  points  using  third-order 
CWENO  reconstruction,  (3.1),  and  Shu’s  SSP  Runge-Kutta  solver,  (3.19). 


algorithms.  In  both  cases,  the  time  scale,  At  is  determined  dynamically  with  CFL=0.475, 


At  = 


_ 0.475 _ 

sj (maxfc  |afc(u)|/Aa:)2  +  (maxfc  \bk(u)\/Az)2 


(5.4) 


where  {ak(u)}k  and  {5fc(w)}fc  represent  the  eigenvalues  of  the  Jacobian  matrices  of  f[u)  and  g(u)  respectively. 


5.1  Transverse  Kelvin— Helmholtz  Instability 

The  Kelvin-Helmholtz  instability  arises  when  two  superposed  fluids  flow  one  over  the  other  with  a  relative 
velocity.  It  models,  for  example,  the  important  mechanism  for  the  momentum  transfer  at  the  Earth’s  magne¬ 
topause  boundary,  which  separates  the  solar  wind  from  the  Earth’s  magnetosphere,  [14].  We  approximate  the 
solution  of  he  two  dimensional  periodic  and  convective  models  with  transverse  magnetic  field  configuration  with 
the  semi-discrete  scheme  of  Kurganov  and  Tadmor,  (2.24)-(2.26),  implemented  along  with  the  reconstruction 
of  §3.1  and  the  SSP  Runge-Kutta  solver  (3.20). 

In  both  cases,  the  governing  equations  are  (5.1)-(5.3)  are  subject  to  initial  conditions 

(p,vx,vy,vz,Bx,By,Bz,p)T  =  (1.0,  vxo  +  vxo,  0, 0, 0, 1.0, 0, 0.5)T,  (5.5) 


where 


VxO 


VxO 


and 


f  -u0sin(^)T^, 

l  o, 


if  —  |  <x  <  | 
otherwise, 


(5.6) 

(5.7) 


with  vo  =  2,  Vo  =  —0.008,  A  =  57t  and  a  =  1.  Also,  the  grids  are  stretched  in  the  2— direction  with  a  Roberts 
transformation,  [14], 


H  sinh  (tz/2H) 
sinh  (r/2) 


t  =  6, 


(5.8) 
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which  renders  a  denser  grid  near  z  =  0  where  the  effect  of  the  small  initial  perturbation  vxq  is  more  noticeable, 
and  a  coarser  grid  near  2  =  ±H ,  where  little  action  takes  place. 

In  the  periodic  case,  the  computational  domain  is  [—  4]  x  [0,  H],  with  L  =  and  H  =  10.  The  free  outflow 

condition  is  applied  at  the  top  boundary,  z  =  H,  and  the  bottom  boundary  values  are  recovered  by  symmetry, 
since  p,  p,  and  By  are  symmetric,  and  vx  and  vz  are  antisymmetric  under  the  transformation  (x,  z)  — >  (— x,  —z). 
In  figure  (5.1),  we  present  solutions  at  t  =  144,  with  72  x  45  and  96  x  60  grid  points. 

The  resolution  and  accuracy  of  our  results  are  comparable  to  those  reported  in  [4]  obtained,  over  two  grids  of 
sizes  96  x  96  and  a  192  x  192  respectively,  with  the  second-order  staggered  central  scheme  of  Jiang  and  Tadmor; 
the  gas  kinetic  scheme  of  Tang  and  Xu,  [29],  using  a  200  x  200  uniform  grid,  and  to  those  obtained  with  Jiang 
and  Wu’s  fifth-order  WENO  scheme,  [14],  over  two  grids  of  sizes  48  x  30  and  96  x  60,  respectively.  The  higher 
computational  cost  of  the  third-order  reconstruction  as  compared  to  that  of  the  fully-discrete  second-order 
scheme  we  used  in  [4]  is  compensated  by  the  fewer  number  of  grid  points  needed  to  resolve  accurately  the  steep 
gradients  that  characterize  the  solution. 


Figure  5.1:  Results  of  transverse  Kelvin-Helmholtz  instability  with  periodic  x— boundary  conditions  computed 
with  third  order  scheme.  Left  column  uses  72  x  45  points  and  right  column  uses  96  x  60  grids  respectively. 
There  are  20  contours  for  density  and  pressure.  Red-high  value,  blue-low  value.  Density  ranges  from  0.79  to 
1.2,  pressure  range  from  0.32  to  0.71  and  maximum  value  for  the  velocity  is  1.25. 

In  the  convective  setup,  the  initial  conditions  and  perturbation  are  the  same  as  in  the  periodic  setup,  (5.5)- 
(5.7).  In  this  case,  the  free  outflow  condition  is  applied  also  in  the  x-direction  over  the  computational  domain 
[—  -j,  -j]  x  [0,  H]\  where  H  =  20  and  L  =  557 r,  with  L  »  A-so  chosen  to  allow  the  excitation  to  convect  freely 
without  disturbing  the  x— boundaries.  The  values  of  the  bottom  boundary  of  the  computational  domain,  2  =  0, 
are  recovered  by  symmetry  as  in  the  periodic  configuration. 

Figures  5.2  and  5.3  display  the  solution  over  the  region  [—50,50]  x  [—20,20]  computed  with  792  x  92  grid 
points  over  computational  domain  at  t  =  120  and  t  =  145  respectively.  The  values  of  the  bottom  half  of  the 
solution  domain  are  recovered  using  the  same  symmetry  conditions  as  used  to  reconstruct  the  bottom  boundary. 
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Figure  5.2:  Solution  of  convective  Kelvin-Helmholtz  instability  at  t=120,  computed  with  third  order  scheme  on 
a  792  x  72  grid.  The  density  ranges  from  0.63  to  1.3,  and  the  pressure  ranges  from  0.20  to  0.85,  the  maximum 
value  for  the  velocity  is  1.54. 
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Figure  5.3:  Solution  of  convective  Kelvin-Helmholtz  instability  at  t=145,  computed  with  third  order  scheme  on 
a  792  x  72  grid.  The  Density  ranges  from  0.43  to  1.3,  pressure  ranges  from  0.10  to  0.86  and  maximum  value  for 
the  velocity  is  1.94. 


Figure  5.4:  Time  evolution  of  the  total  transverse  kinetic  energy,  log(^  J  pv^dxdz),  integrated  over  [-L/2,  L/ 2]  x 
[—H,H],  for  both  periodic  and  convective  Kelvin-Helmholtz  instability.  The  results  for  the  periodic  case  with 
96  x  96  and  192  x  192  grid  points  are  represented  by  a  dashed  and  a  doted  curve  respectively.  The  convective 
configuration  is  represented  by  a  solid  line. 
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5.2  Orszag-Tang  MHD  Turbulence  Problem 

This  model  considers  the  evolution  of  a  compressible  Orszag-Tang  vortex  system.  The  evolution  of  the  vortex 
system  involves  the  interaction  between  several  shock  waves  traveling  at  various  speed  regimes  [14,  32],  which 
makes  the  problem  useful  to  validate  the  robustness  of  numerical  schemes. 

The  initial  data  is  given  by 

p(x,  z,  0)  =  y2,  vx(x,  z,  0)  =  —  sinz,  vz(x,  z,0)  =  sin  a;, 
p(x,  z,  0)  =  7,  Bx{x,  z,  0)  =  —  sin  z,  Bz(x,  z,  0)  =  sin  2x 

where  7  =  5/3.  With  this  data,  the  RMS  values  of  the  velocity  and  magnetic  fields  are  both  1;  the  initial 
average  Mach  number  is  1,  and  the  average  plasma  beta  is  10/3. 

We  solve  the  problem  in  [0,  27r]  x  [0,  27t],  with  periodic  boundary  conditions  in  both  x-  and  z-directions.  For 
this  problem  we  implement  the  fourth-order  reconstruction  of  §3.2  and  the  SSP  Runge-Kutta  fourth-order  ODE 
solver  (3.20)  using  a  uniform  grid  with  192  x  192  points. 

In  this  problem,  the  numerical  scheme  fails  to  satisfy  the  divergence  free  constraint  of  the  magnetic  field, 
V  •  B  =  0.  In  order  to  guarantee  this  condition  and  avoid  numerical  instability,  we  project  the  updated  magnetic 
field  B  into  its  divergence  free  component  at  the  end  of  every  time  iteration  by  applying  the  correction  (3.21)- 
(3.22). 


Figure  5.5:  Fourth-order  solution  of  Orszag-Tang  MHD  turbulence  problem  with  a  192  x  192  uniform  grid  at 
t  =  0.5.  There  are  20  contours  for  density  and  pressure.  Red  -high  value,  blue  -low  value.  Density  range  from 
2.1  to  5.9,  pressure  range  from  1.0  to  5.8.  The  maximum  values  of  |u|  and  \B\  are  1.6  and  1.6  respectively. 

Figures  5. 5-5. 7  display  the  solution  of  the  Orszag-Tang  vortex  system  at  t  =  0.5,  t  =  2,  and  t  =  3  respectively. 
These  simulations  were  performed  with  192  x  192  grid  points  using  the  fourth-order,  genuinely  multidimensional, 
CWENO  reconstruction  of  §3.2  SSP  Runge-Kutta  solver  (3.18)-(3.20).  These  results  demonstrate  the  ability 
of  higher-order  central  schemes  to  resolve  the  shocks  that  the  vortex  system  develops  while  maintaining  the 
simplicity  and  ease  of  implementation  typical  of  this  black-box  type  of  finite-difference  schemes.  We  also  note 
that  similar  results  were  obtained  with  the  dimension-by-dimension  third-order  reconstruction  and  the  Runge- 
Kutta  solver  (3.19)  using  a  288  x  288  mesh.  The  inproved  resolution  of  the  fourth-order  scheme  allows  us  to 
compute  accurate  approximations  using  coarser  grids  than  those  required  by  lower  order  schemes,  thus  reducing 
the  computational  cost  of  this  type  of  simulations;  an  advantage  that  is  more  evident  in  higher  space  dimensions. 
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Figure  5.6:  Second-order  solution  of  Orszag-Tang  MHD  turbulence  problem  with  a  384  x  384  uniform  grid  at 
t  =  2.  There  are  20  contours  for  density  and  pressure.  Red-high  value,  blue-low  value.  Density  ranges  from 
0.62  to  6.3,  pressure  ranges  from  0.14  to  7.0.  The  maximum  values  of  |n|  and  \B\  are  1.6  and  2.8  respectively. 


Figure  5.7:  Fourth-order  solution  of  Orszag-Tang  MHD  turbulence  problem  with  a  192  x  192  uniform  grid  at 
t  =  2.  There  are  20  contours  for  density  and  pressure.  Red  -high  value,  blue  -low  value.  Density  range  from 
0.62  to  6.3,  pressure  range  from  0.14  to  7.0.  The  maximum  values  of  |u|  and  \B\  are  1.6  and  2.8  respectively. 
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Figure  5.8:  Pressure  distribution  along  the  line  z  =  0.6257T  at  t  =  3.  Dashed  line  corresponds  to  third-order 
approximation  with  288  x  288  grid  points  and  dotted-line  corresponds  to  fourth-order  results  computed  with 
192  x  192  grid  points. 
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